Morphological multiparameter filtration and persistent homology in mitochondrial image analysis

The complexity of branching and curvilinear morphology of a complete mitochondrial network within each cell is challenging to analyze and quantify. To address this challenge, we developed an image analysis technique using persistent homology with a multiparameter filtration framework, combining image processing techniques in mathematical morphology. We show that such filtrations contain both topological and geometric information about complex cellular organelle structures, which allows a software program to extract meaningful features. Using this information, we also develop a connectivity index that describes the morphology of the branching patterns. As proof of concept, we utilize this approach to study how mitochondrial networks are altered by genetic changes in the Optineurin gene. Mutations in the autophagy gene Optineurin (OPTN) are associated with primary open-angle glaucoma (POAG), amyotrophic lateral sclerosis (ALS), and Paget’s disease of the bone, but the pathophysiological mechanism is unclear. We utilized the proposed mathematical morphology-based multiparameter filtration and persistent homology approach to analyze and quantitatively compare how changes in the OPTN gene alter mitochondrial structures from their normal interconnected, tubular morphology into scattered, fragmented pieces.


Introduction
Recent advances in cellular biology have revealed that the majority of intracellular organelles are not distinct units with regular geometric shapes.They appear as complex, irregular threedimensional (3-D) structures with complex branching patterns or "holes" within structures devoid of organelle material.For example, mitochondria are more than individual ovoid-like structures but are found within an extensive network of curvilinear and branching tubular structures.Smaller mitochondria pieces are continuously separated from the main mitochondrial network in a process called fission when senescent or damaged in order to be degraded.On the other hand, two distinct mitochondria can combine and merge their inner and outer mitochondrial membranes in a process called fusion.When performing quantitative analysis of the shape and size of these organelles, it is difficult to study morphological changes in the overall geometry of the irregular mitochondrial network between cells.This challenge is exacerbated by differences in the shape of cells, which in turn influences organelle geometry.Current analytic methodologies are limited in their ability to take into account the overall irregular shapes for quantitation and to compare differences between cells.Therefore, the purpose of our work is to describe our novel mathematical representation of complex biological organelle shapes using persistent homology based on the multiparameter filtration [1,2] and demonstrate the utility of our analysis methodology in analyzing morphological differences in intracellular organelles such as mitochondria.
In order to study cells and quantify organelle size and density at the histological level, cellular tissue or a monolayer of cells is often immunostained with a fluorescent dye, mounted on a slide, and imaged with a fluorescence microscope.The advent of the confocal microscope allows cells to be imaged in a thin optical "section" through various parts of the cells.Multiple optical histological sections can be obtained in exquisite detail.Since fluorescence intensity is proportional to the amount of organelle protein or material stained, the fluorescence signal is quantitative and can be analyzed by software for comparison between different samples.
Current software analysis of fluorescence staining of organelles is imprecise due to major limitations.First, to analyze fluorescence signals resulting from stained organelles, current approaches involve a "thresholding" step in which the user or the software determines the fluorescence signal intensity level that separates background noise from actual biologically relevant data.The software then uses the threshold level to calculate the total fluorescence signal and intensity for quantitation.This thresholding process is often performed manually by the researcher, and so can be quite arbitrary, subjective, and highly variable.Different threshold settings result in variability in the quantitation and subsequent comparative image analysis.This problem is exacerbated by high background staining or low fluorescence signal intensity.Secondly, it is important for the software to be able to determine cellular boundaries.Often, to allow relative comparison of organelle abundance among cells of different sizes, the fluorescence intensity of the organelle would need to be normalized by dividing by the cellular area.This process often requires human users to manually "draw" these cellular boundaries for the software, which introduces another subjective variability.Lastly, due to inherent variation in immunostaining experiments, one set of samples may be stained more intensely than other sets of experiments performed on a different day.What may appear to be noise in one set of experiments might appear as an actual fluorescence signal in another experiment.This often exacerbates the other two limitations mentioned above.
To address these limitations, we propose an approach to mathematically model a fluorescently stained organelle using both persistent homology (PH) and a multiparameter filtration (MF) method, aiming to reduce experimental and data analysis variability.PH, MF, and the induced multiparameter persistent homology (MPH) belong to the burgeoning field of Topological Data Analysis (TDA) (see, for instance, [2]).Recently, the PH framework has demonstrated its potential in microscopy image analysis, including applications in pathological bone [3], cancer-immune microenvironment [4], and the deep learning framework in medical image segmentation [5,6].On the other hand, the concept of MPH, representing a more general version of PH, has been well-established [7,8], with studies applying MPH in diverse areas, such as biomolecular data [9][10][11] and biomedical/biological image analysis [12].Notably, research incorporating MPH into the analysis of microscopy images remains limited and rare.In this work, we propose a hybrid framework that considers PH features on MF, namely, the MF-PH framework, incorporating operations in mathematical morphology to encode local geometric and topological information within digital mitochondria images.Compared to conventional PH based on sub-level set filtrations, this approach offers more

Cell cultures
Primary mouse embryonic fibroblasts (MEFs) were purified from wild-type (WT), optineurin (OPTN), and OPTN knock-out (KO) transgenic mice previously generated in our laboratory [20,21].MEF cells from WT-OPTN and KO-OPTN mice were plated on sterile glass coverslips in 24 well plates.They were cultured at 37˚C in a humidified 5% CO 2 incubator in media consisting of DMEM that was supplemented with 10% FBS, and 10% penicillin-streptomycin and GlutaMAX (all culture media reagents were from Gibco/ThermoFisher).The cells were passaged between 75-85% confluence to maintain the cell line.Cultures were allowed to sit overnight after plating before any manipulation of the cultural environment was performed.Cells were fixed in 4% paraformaldehyde.

Immunofluorescence staining
Following fixation, cells were washed with PBS, and incubated with a primary antibody of Tom20 (Santa Cruz Biotechnology) for mitochondria in 5% donkey serum and 5% Triton X-100 in PBS overnight.After several PBS washes, appropriate secondary antibodies conjugated to Alexa Fluor-488 (1:500) (Invitrogen) were added to the coverslips for 1h before mounting on coverslips with DAPI Vectashield mounting medium (Vector Labs).

Confocal fluorescence microscopy
In the experiment of the work, two inhibitors for the mitochondrial cell are used, oligomycin and antimycin A, which impair mitochondrial respiration to induce mitochondrial damage and mitophagy.Especially, treatment with oligomycin and antimycin A for longer amounts of time should result in greater mitochondrial damage, subsequently leading to more fragmented and punctate mitochondrial networks.In this work, we consider two classes: WT and KO.Each class contains five images.In addition, we consider each class with oligomycin and antimycin A for 2 hours and 4 hours, respectively.The total images we have are 30 (5 × 3 × 2) images.Fig 1 demonstrates samples of confocal images.All cells were imaged with the same microscopy settings.However, fluorescence signals in KO cells appear lighter due to fragmentation of the mitochondrial network and decreased abundance of mitochondria proteins, particularly the Tom20 protein that was used as a marker for immunostaining.Mitochondrial network fragmentation is due to the genetic knockout of the optineurin gene.
Cells were imaged with a Nikon C2 Plus confocal microscope using a 63x NA oil immersion objective.A 488 nm laser was used to excite AlexaFluor 488.Slides were imaged at a resolution of 2048 pixels by 2048 pixels with the pinhole set to 1.2 AU.An exposure time of 1.2 μs per pixel was used, and the average of two frame acquisitions was used as the raw microscope image.Regions of interest were selected and a series of z-sections were taken from the top to the cell bottom with intervals between sections set to 2 μm.

Ethical approval
All mouse experiments were performed at Duke University after having received approval from the Institutional Animal Care & Use Committee (IACUC) and were conducted according to Association for Research in Vision and Ophthalmology guidelines.

Methods
The goal of this work is to use multiparameter filtration-based one-persistent homology to analyze confocal images of mitochondrial networks.This section has two main objectives: first, to review the essential mathematical details needed to introduce the proposed methods, and second, to present the mathematical notations that will be used throughout the article.This section begins with a review of fundamental topics related to digital images, including the use of binary and grayscale images and the mathematical modeling of complex shapes observed in mitochondrial morphology.It then covers the concept of Betti numbers and how they are assigned to a given binary image.Following this, we introduce and apply the concept of filtration, and discuss persistent homology.An overview of notations is provided in Table 1.For more mathematical background, we refer readers to these publications [1,7,[22][23][24].

Digital images
Digital images consist of discrete grid points called pixels.Each pixel is a square and has a coordinate.The set of 2-dimensional integer lattice Z 2 is a common mathematical object to model the coordinate of pixels.Since digital images have fixed width and height, we use P which is a subset of Z 2 (P � Z 2 ) to denote the collection of pixels.For instance, Fig 2 shows a 12-by-12 binary image, and thus, P is the set of those 12-by-12 grids (pixels).
A binary image is then a function f : P � Z 2 !f0; 1g, where 0 (1) usually represents the black (white, respectively) color.In this work, our focus will be the shapes formed by black pixels, i.e., the collection of pixels whose value is 0, formally denoted by K(f) ≔ {(x, y) 2 P | f(x, y) = 0}.Interpreting the black pixels as black cubes within the R 2 plane, K(f) can be recognized as a subspace of R 2 , commonly referred to as a cubical complex [25].Reviewing Opening operation on f with the structural element S i .
The number of isolated black regions in f.
The number of white regions enclosed by black regions in f.
Filtration: a collection of sets X i with the property Sublevel set filtration of a grayscale image g.
Opening filtration of a binary image f.
The the left side illustrates the visualization of a binary image, while the right side depicts the corresponding function values of the binary image f.Moreover, in this example, K(f) represents the black regions (shape of "1" and "0").
Grayscale images are the extension of binary images.Confocal images of the mitochondrial network in this study are grayscale images.Formally, an 8-bit grayscale image is a function g : P � Z 2 !f0; 1; 2; � � � ; 255g.The main difference between a grayscale and a binary image is that the pixel values of a grayscale image range from 0 to 255 while those of a binary image range are either 0 or 1.  Converting a grayscale image to a binary image has been described in computer vision [26,27].In this work, we consider the most common way of doing this: a global thresholding method, which means that pixel values below a certain threshold are set to 0 and 1 otherwise.In other words, pixel values below a certain threshold are set to black and white otherwise.Mathematically, it can be described as follows: Given a grayscale image g and a threshold value t, the thresholded binary image of g at t is defined as Consequently, given an 8-bit grayscale image g, global thresholding produces a total of 256 binary images: show the binary image at different threshold values.Selecting one threshold is a challenging task and is often subjective (see e.g.[28][29][30].On the other hand, our proposed method takes all threshold values into consideration.At this point, we have seen binary, grayscale images, and their interaction via the global thresholding method in (1).Digital image processing concerns the manipulation of binary or grayscale images.Among many tools in the field, we consider a classic topic Mathematical Morphology that offers a wide range of methods to analyze images [31][32][33][34][35][36][37][38].
Methods in mathematical morphology are often called morphological operations.The building block of morphological operations is the structuring element.Because morphological operations are concerned with interactions among pixels, in particular local information around given pixels, structuring elements are the medium to acquire such information.We consider the opening operation associated with the square structuring element in this work.We denote it by O S i ðf Þ, where f is a given binary image, and S i is the square structuring element of size (i + 1) × (i + 1).Specifically, S i can be depicted below where the white bullet � denotes the original point (0, 0) in Z 2 .One may think of the opening operation as a transformation that takes the given binary image and outputs the processed binary image.In other words, given a binary image f and a structural element S, applying the opening operation O S to f results in another binary image, denoted by O S (f).Note that when S = S 0 , the opening operation O S is identity, i.e., O S 0 ðf Þ ¼ f , and that O S i ðf Þ is also a binary image.Its formal mathematical definition can be found in these publications [1,[31][32][33].Opening operation is commonly used to reduce noises in the binary image.By its design, an opening operation removes isolated white regions whose size is less than the structural element S. For instance, Fig 4(a) shows a binary image with various isolated white regions.We see that in

Betti numbers
The Betti numbers constitute a classic subject in the field of algebraic topology in mathematics [39][40][41][42].They play a crucial role in quantifying objects based on their connectivity, as well as in characterizing loop and void structures.They have been used in digital image processing [43][44][45][46][47]. Recently, the Betti numbers of one-parameter persistent homology were discovered to exhibit capabilities in identifying and segmenting tasks within microscope image analysis [48].In this subsection, we will explain the meaning of Betti numbers for binary images.
Let f be a binary image.We consider its black region, K(f).Betti numbers of K(f), denoted by β(K(f)), is a pair of numbers (β 0 , β 1 ) where β 0 (β 1 ) is called 0-th (1-st) level Betti numbers.β 0 counts the number of isolated black regions, and β 1 counts the number of white regions that are completely enclosed by the black ones.Take Thus, setting the manual threshold is often subjective, and differences in thresholding will change the analysis outcome.To address this challenge, our proposed MF approach, on the

Filtration
At the end of the previous subsection, we observe that different thresholds or different structural elements would result in different binary images.To address this, we will utilize the Persistent Homology (PH), a recently developed tool in the field of topological data analysis [7,[49][50][51][52][53][54].We first introduce the concept of filtration, which is essential for understanding persistent homology.
Filtration is a sequence of sets where sets satisfy subset relations.Mathematically, let X i be some set with index i.In our context, we think of X i as the black set of the binary image.Consider a family of sets fX i g n i¼0 and we say that The definition means that the set X i grows as the index i increases.In our context, since we consider the black sets of binary images, filtration of black sets means the binary image becomes darker as the index increases.There are two filtrations in this work: sublevel set filtration and opening filtration.
The first filtration is related to the global thresholding method in (1).Given an 8-bit grayscale image g, recall from (1), g t is the thresholded binary image at threshold t.We consider the set of black pixels of g t , i.e., K(g t ).The sublevel set filtration (cf.[53], Sec.VI.1, p. 125) is the following nested subset relations: The idea is that as the threshold value increases, the corresponding binary image becomes darker (in the extreme case, when the threshold value is 255, then the binary image is The second filtration is related to the opening operation.As discussed in the previous subsection, opening filtration with the structural element S i can be thought of as filtering those isolated white regions whose size is less than (i + 1) × (i + 1).Thus, as S i becomes larger, more white regions will be removed by the opening operation; for instance, it can be seen in Fig 4.
We now consider a family of opening operations O S i for i = 0, 1, 2, � � �, 20.As shown in [1,55], this family of operations possesses the subset relations called the opening filtration: where f is a given binary image.The opening filtration of f is diagram (5) and is denoted by

Persistent homology
Here, we introduce the idea of persistent homology [51,53,54] and its representation: persistence diagrams.At this point, we have seen Betti numbers of the binary images.Persistent homology can be thought of as the generalization of the Betti numbers.Formal mathematical developments of persistent homology can be found in [51,53,54].In this subsection, we will use an example to demonstrate persistent homology.Filtration of sets is the essential assumption for persistent homology.Given a filtration fX i g n i¼0 persistent homology not only computes Betti numbers of each set but also tracks the changes of Betti numbers.To see persistent homology in action, we demonstrate it via a concrete example of a grayscale image below.(6) The corresponding threshold decomposition, i.e., all possible binary images of g, are shown in Fig 7.
To determine persistent homology, we calculate Betti numbers at each thresholded binary image K(g t ) and track changes in Betti numbers because of the filtration relation.Note that the Betti numbers for each K(g t ) can be found in Fig 7 and also that there are three changes of Betti numbers, i.e., (1) from K(g 1 ) to K(g 2 ) where Betti numbers change from (1, 0) to (2, 0), (2) from K(g 2 ) to K(g 3 ) where Betti numbers change from (2, 0) to (1, 1), and (3) from K(g 9 ) to K(g 10 ) where Betti numbers change from (1, 1) to (1, 0).Persistent homology keeps track of these changes in Betti numbers.
On the ground level, we will look into how pixels form components/holes, or fill up the grid.In Fig 7, since at K(g 1 ) there is a newly formed black column, we say that "a β 0 feature is born at K(g 1 )".Similarly, at K(g 2 ), there is another newly formed black column, and another β 0 feature is born at K(g 2 ).At K(g 3 ), two changes occur.One is that the β 0 feature that is born at K (g 2 ) joins with the other β 0 feature.Thus, we say that this generator dies at K(g 3 ).For this generator, we record its birth and death time as a pair of numbers (2,3), where the first (second) component is the birth (death) value.We also observe that at K(g 3 ), β 1 changes from 0 to 1, and thus, a β 1 generator is born at K(g 3 ).From K(g 4 ) to K(g 9 ), black sets do not change.Finally, at K(g 10 ), the entire image becomes black and thus the Betti number is (1, 0).A β 1 feature dies at K(g 10 ), and its birth-death pair would be (3,10).Note that the first β 0 generator persists in this process, so by convention, we assign its death value by 1 (cf.[53], Sec.VII.1, p. 152), its birth-death pair is (1,1).Collections of such birth-death pairs are called persistence diagrams with respect to the filtration.The persistence diagrams associated with the dimensions of connected components and loops are termed the persistence diagrams in dimensions 0 and 1, respectively.We use the following notation: P 0 ðfKðg t Þg 10 t¼1 Þ ¼ fð1; 1Þ; ð2; 3Þg and P 1 ðfKðg t Þg 10 t¼1 Þ ¼ fð3; 10Þg.To summarize this subsection, given a filtration of sets fX i g n i¼0 , there is a corresponding persistence diagram denoted by PðfX i g n i¼0 Þ.There are several sophisticated software packages to compute persistence diagrams for a given filtration (see a survey article [56] for a list of software).Among them, we use Perseus software package [57,58] to compute persistence diagrams in this work.

Multiparameter filtration construction
At this point, we have seen that given a filtration, one can study persistence diagrams.In particular, such filtration is with respect to one parameter, such as threshold values in sublevel set filtration in (3) or the sizes of structural elements in the opening filtration (5).It is also known as one-parameter filtration, and the corresponding persistence diagrams are called one-parameter persistence.Recently, there have been studies to generalize this concept to multiparameter filtration and multiparameter persistence [22][23][24][59][60][61].The proposed model will utilize two-parameter filtration or bifiltration.In what follows, we will introduce the proposed model.The core idea of the proposed model is to combine sublevel set filtration in (3) and opening filtration in (5).
We are now ready to present our proposed model.Given a grayscale image g, the proposed model combines sublevel set filtration in (4) and opening filtration in (5) to form a bifiltration.The rigorous consideration and its mathematical proof can be found in [1].Intuitively, one may start with the sublevel set filtration Kðg i 1 Þ � Kðg i 2 Þ as long as i 1 � i 2 .Observe that since for fixed i, g i is a binary image, one may apply opening operations to g i to obtain the opening filtration KðO S j 1 ðg i ÞÞ � KðO S j 2 ðg i ÞÞ for any j 1 � j 2 .Therefore, for each binary image g i , there is a corresponding opening filtration; moreover, each g i satisfies the sublevel set filtration.Thus, we obtain the following bifiltration: The Eq (10) is our proposed persistent homology model based on multiparameter filtrations.Given a grayscale image g, we then consider the corresponding bifiltration.Figs 8 and 9 visualize the bifiltration in (10).
As discussed in the previous sections, given a filtration of sets, one may track changes in Betti numbers of these multisets in the form of persistence diagrams.Similarly, one would desire to possess analogous "persistence diagrams" for a multiparameter filtration.Unfortunately, the notion of persistence diagrams for a multiparameter filtration does not exist [8,62], and hence, it is challenging to utilize its full power in practice.There have been works to devise ways for practical use.For example, the software RIVET [63] offers a framework for constructing a bifiltration for a given point cloud data and computes a 1-dimensional persistence diagram along any path in the bifiltration [12] built upon RIVET to develop a multiparameter landscape which can be seen as extracting information from the bifiltration.In our study, we aim to extract features from the bifiltration as well.The main difference is that our approach is specifically for image data (cubical complexes) as opposed to point cloud data (simplicial complexes).In addition, the ways of constructing bifiltration as we describe in this article are by morphological operations.In the next section, we will describe our proposed features from the proposed bifiltration.We summarize this section by stating the proposed pipeline: Beyond bifiltrations, the proposed MF construction based on morphological operations can be generalized to MF with more parameters.For instance, suppose S 1 � S 2 � � � � � S n and T 1 � T 2 � � � � � T n are two increasing sequences of structuring elements as required in [1,55] Actually, by leveraging different morphological operations (e.g., erosion, dilation, closing, etc. [32]), a multifiltration of images in any number of parameters can be produced [1].However, to demonstrate multifiltration in a more geometrically intuitive and interpretable way, we focus on investigating the bifiltration described in (10) and compute the persistent homology on it.
The extraction of the PH information from the proposed bifiltration proceeds as follows.Note that each horizontal relation in ( 10) is a sublevel set filtration; each vertical relation in ( 10) is an opening filtration.Each horizontal contains topological information: how the grayscale image is formed as the threshold increases.Each vertical contains geometric information: the size of holes at each binary image.It is also possible to follow a path that is neither horizontal nor vertical, which could be an interesting topic for future research.In this study, however, we will focus on vertical and horizontal paths.More precisely, we consider P H;i 1 is the 1-st persistence diagram of the horizontal path for fixed i; P V;t 1 is the 1-st persistence diagram of the vertical path for fixed t.

Proposed features
As shown in the previous section, diagram (10) offers rich information about both the geometric and topological properties of the underlying images.In this section, we will describe the features from the bifiltration in (10) and apply them to the mitochondrial images.Each horizontal or vertical filtration contains different information about mitochondria.There are three main features that we will use in this work: (A) normalized Betti numbers curve, (B) size distribution of the Betti numbers, and (C) the connectivity index that we present.We will describe details of each feature in this section.These three features are summarized in Table 2.

Normalized Betti curve
Betti curves have been used in several studies and shown to be effective shape descriptors [64][65][66][67].Normalized Betti curves, on the other hand, are more recently introduced in the work by [65].The essential idea is to calculate β 1 for each thresholded binary image at t and normalize it by the total number of points in a persistence diagram.Compared to the classic Betti number curve, [65] showed that it possesses desired theoretic properties over the classic one.We refer interested readers to [65] for more details.Building upon this work, we are using normalized Betti curves in the context of bifiltration for our MF-PH analysis technique.
To illustrate our approach, let the opening parameter i be fixed and consider the horizontal persistence diagram in (15), P H;i 1 .The normalized Betti curve is defined for the i-th row of (10) as where denotes the number of elements in a set or a multiset.By its design, we observe that the values of r H;i 1 are between 0 and 1.Therefore, for any different raw image, this feature is both scale and size-independent.No preprocessing on the raw images is required in order to compare these curves.

Size distribution of Betti numbers
Betti numbers are invariant under scaling, so objects with different sizes would count equally toward Betti numbers [39][40][41][42].Because of our proposed bifiltration (10), such geometric information may be captured from the vertical filtration.This idea is similar to the so-called granumonetry [32,[68][69][70] in material science.
Our rationale is that since mitochondrial networks in KO cells are fragmented, these are expected to be smaller in size.In contrast, WT cells have mitochondrial networks that are not as fragmented, so consist of mitochondria of larger sizes.To quantitatively compare differences in sizes, this information is represented by the vertical direction in (10).We designed this feature based on the size of Betti numbers, and thus we will call it the size distribution of Betti numbers.For each t, we consider the vertical direction in (10) and the corresponding persistence diagram in (16) P V;t 1 .As studied in [1], the multiset of death values D t ¼ fd j ð0; dÞ 2 indicates the sizes of 1-st dimensional holes in the binary image X t,0 .Therefore, the size distribution of the Betti numbers along t-th column of diagram ( 10) is defined as Since our opening filtration is from 1 to 20, the number of bins in ϕ(D t ) is 20.For instance, at the 1-st bin of ϕ(D t ), it represents the percentage of Betti numbers of X t,0 that are of size 1 × 1; at the 2-nd bin of ϕ(D t ), it represents the percentage of Betti numbers of X t,0 that are of size 2 × 2; at the 20-th bin of ϕ(D t ), it represents the percentage of Betti numbers of X t,0 that are of size 20 × 20.

Connectivity index
Mitochondrial fragmentation and loss of connectivity are crucial elements in the analysis of how mitochondrial morphology is altered in response to genetic, biological, pharmacological, or environmental changes.Therefore, we propose a feature called connectivity index to assess the connectivity of a mitochondrial network.By identifying D t as a sub-multiset of P V;t 1 via the mapping d 7 !(0, d), we design our connectivity index as Note that by definition, C t is calculated along the t-th column filtration in (10).Furthermore, the porous structures with persistence intervals of the form (b, 1) are excluded from the connectivity index computation.By considering the set size of D t and the lifespans d − b of the persistence intervals, the proposed connectivity index offers a description of the fragmentation and loss of connectivity of a given image.By definition, the values of C t range from 0 to 1 for any fixed t.If C t is closer to 1, the input image is more connected; on the other hand, if C t is closer to 0, the input image is more fragmented.As mentioned in (18), D t indicates the sizes of Betti numbers in the binary image X (t,0) .If D t ¼ P V;t 1 , then the size of 1-dimensional holes (or porous structures) in the binary image X t,0 can be determined exactly.For otherwise, P V;t 1 n D t 6 ¼ ; means there are additional holes that are coming off the original holes.In other words, the more features in the multiset P V;t 1 n D t , the more connected the original binary image X (t,0) .Fig 10 illustrates examples of three binary images to interpret the geometric insight of the connectivity index on digital images.The images located in the front rows have more sparse and fragmented (white) structures.As Betti-1 information, small porous structures as white pixels disappear as the black pixels increase in the opening filtration and correspond to smaller death values in the persistence diagram P V;t 1 .Furthermore, Betti-1 structures encoded as members in the set D t have the birth value 0, representing the original structure in the image.For instance, for the first binary image in Fig 10, numerous porous structures exist in the original image, belonging to D t , leading to a smaller connectivity index (0.43).On the contrary, the third binary image is more robust and contains fewer singleton pore structures.Compared to the first and the second images, the cardinality of D t is much smaller, leading to a higher connectivity index (0.9).

Computational setup
We chose Perseus software to compute persistence diagrams because it offers executable files for all major operating systems (Windows, Linux, and Mac) and can be adapted to any coding environment [57].Updates of the proposed framework in other packages (e.g., Gudhi [71]) can be found in our GitHub repository (see the "Code availability" section).
Using Perseus, computing a single persistent homology for the persistence diagram P H;i 1 or P V;i 1 for the proposed curve (r H;0 1 ðtÞ, ϕ(D t ), or C t ) can be done in 10 to 20 seconds.The computation is based on a Windows 11 x64 environment with an AMD Ryzen 7 7840HS CPU and 16.0 GB of RAM.

Results
This section presents the results for the proposed features on mitochondria images of types KO and WT: Normalized Betti Curve, Size Distribution of Betti Numbers, and Connectivity Index.The following paragraphs provide a detailed analysis of each feature. 1 , r H;3 1 , and r H;5 1 , respectively.We observe that in (b), the shapes of these curves are similar among different sample images.Similar behaviors can also be observed in (e) and (h).The similarity in these curves suggests that all of these KO OPTN cells exhibit mitochondrial morphologies that are similar.

Normalized Betti curve
Lastly, we compare the curves between WT and KO types.To facilitate the comparison, we consider curves in Fig 11(c), 11(f) and 11(i) in which we generated curves representing averages for all the WT and KO cell samples.For instance, in (c), the blue curve represents the average curve over those 5 sample curves in Fig 11(a) while the red curve represents the average curve over those 5 curves in Fig 11(b).Similarly, the blue curves in (f) and (i) represent the average curve over those curves in (d) and (g), respectively.We observe that the difference in (c) is not obvious.However, in (f) and (i), the average curves for WT and KO are quite different.(c), (f), and (i) represent morphological comparisons across all gray-level thresholding.Therefore, our findings are consistent with what we see visually in these microscope images.This may suggest that normalized Betti curves are good descriptors for these two types of cells.

Size distribution of Betti numbers
We apply the size distribution of Betti numbers in (19) to quantify WT and KO.ϕ(D 128 ), for WT and KO type, respectively.Each curve represents a single-cell sample; there are 5 different cells used per genotype.We observe that in Fig 12(a), among images in WT, the distributions seem to have more variability.For instance, the green curve in (a) shows that about 10% of Betti numbers are of size (1 + 1) × (1 + 1) (the first bin) and about 5% of Betti numbers are of size (7 + 1) × (7 + 1) (the 7-th bin); the red curve in (a) shows that about 50% of Betti numbers are of size (1 + 1) × (1 + 1) and none of the Betti numbers are beyond size (6 + 1) × (6 + 1).This behavior of high variability was also seen in the case of the normalized Betti curve.
On the other hand, in Fig 12(b), we observe that the distributions are very similar among the 5 sample images representing different cells for the KO genotype.For instance, curves in (b) show that more than 50% of Betti curves are of size (1 + 1) × (1 + 1), and almost none of them are of size (5 + 1) × (5 + 1) (zero after the 5-th bin).This reflects our expectation that mitochondrial networks in KO cells are more fragmented.
To further quantify the difference between WT and KO type, we calculated averages over curves in Fig 12(a) and 12(b), respectively.The average curve plots are shown in Fig 12(c) and 12(d), respectively.We observe that for the WT type, mitochondrial sizes of 1, 2, and 3 are about 40%, 20%, and 15%.On the other hand, for the KO type, they are about 70%, 15%, and 5%, which shows that a larger percentage of KO mitochondria are smaller and thus more fragmented than WT mitochondria.The results suggest that WT and KO have different size profiles.Furthermore, we also calculate the total rate of changes  indices C t for WT have more variability.Some cells are more connected than others.Such behavior was also observed in the case of the normalized Betti curve and the size distribution.On the other hand, the connectivity indices for KO cells have small variability because KO mitochondria networks in all the cell samples are quite fragmented.

Statistical analysis
To investigate the feature behaviors of normalized Betti curve, size distribution of Betti numbers, and connectivity index on the sampled datasets shown in Figs 11-13, we performed statistical tests.Specifically, we consider the average values of the Normalized Betti curve and connectivity index curves over the threshold interval from 100 to 150 as the statistical targets for the hypothesis test (cf.Fig 13).On the other hand, based on the size distribution curves along with the sizes (the x-axis) of the structuring elements depicted in Fig 12, we consider the average values of the sampled size distribution curves over the size interval from 5 to 15.Using these average values as the target statistics in the hypothesis tests, permutation tests are performed and the results are depicted in Table 3.The permutation hypothesis tests are performed using the Python package mlxtend [72], with the significance level α set to 0.05.Demonstrations of the permutation hypothesis tests are also provided in our GitHub repository (see the "Code availability" section).
According to the results presented in Table 3, the proposed normalized Betti curves, size distribution of Betti numbers, and connectivity index methods demonstrate their potential in analyzing the differences and mutations between WT and KO mitochondrial structures.Based on the p-values from the experiments (0.009 ± 0.003, 0.010 ± 0.003, and 0.009 ± 0.003), curves r H; 3  1 , r H;5 1 , and C t achieve the most statistical significance for separating the WT and KO classes, showing that he WT and KO classes of images can be effectively separated based on these features.Furthermore, although the experiment of size distribution ϕ(D 128 )(s) yields a higher pvalue and standard deviation (0.025 ± 0.005) compared to curves r H; 3  1 , r H;5 1 , and C t , the result still rejects the null hypothesis that the WT and KO classes share the same mean of the size distributions in the range [5,15].Finally, we observe that the hypothesis test of the curve r H;0 1 doesn't show the statistical significance of the separation of WT and KO classes since the pvalue (0.427 ± 0.017) is larger than the significance level α = 0.05.According to the definition (cf.( 15)), the r H;0 1 curves are produced by computing the conventional persistent homology

Discussion and conclusion
In this work, we propose a novel mathematical framework combining persistent homology, mathematical morphology, and multiparameter filtration to study mitochondrial structural differences.The framework is flexible and can be applied to other types of cells.We also demonstrate ways to extract information from the bifiltration by developing three features for quantitating structural differences-normalized Betti curves, the size distribution of Betti numbers, and the connectivity index.They are independent of the size of each cell and robust to noise.Especially, with a main focus on demonstrating how the mathematical morphology operations integrate into the PH and MF framework, such as opening and thresholding operations, we consider the proposed three features that directly connect to the quantitative structural information within the mitochondria networks.In future investigations, a central focus will revolve around advancing integrated machine learning model with the proficiency to autonomously determine parameters within the proposed features.This includes the adept selection of the sizes and shapes of the structuring elements.Concurrently, an additional avenue of exploration involves the integration of intrinsic features, notably the incorporation of persistence statistics [46], persistence images [47], and persistence landscapes [45].This approach aims to delve deeper into the networks, providing a richer and more nuanced understanding.
Apart from the exploration of one-parameter persistent homology within the 2-dimensional persistence, a future focus of our study is the comprehensive investigation of 2-dimensional persistence.While the RIVET software [63] provides a framework for extracting persistence features and descriptors from a 2-dimensional bilfiltration, it is essential to note that the target filtration is founded on the point cloud within the R n space.This differs from the cubical complex filtration proposed in our work for the image [25].Considering the 2-dimensional algebraic network induced by the bifiltration, a prospective avenue for future research involves a parallel implementation of RIVET's feature extraction.This approach is crucial for aligning with our proposed image-based cubical complex filtration and holds the potential for enhancing the efficiency and relevance of our methodology.Furthermore, [1] proposed a n-parameter filtration on images via cubical complexes (where n > 2) based on alternating morphological operations.It would be interesting to utilize that n-parameter filtration to analyze the mitochondrial networks.
To the best of our knowledge, this work is one of the first ones that utilizes the one-persistent homology based on multiparameter filtrations to quantify mitochondrial networks.Recently, [12] used multiparameter persistence to study images of tumor cells.However, the methods in [12] of treating images and constructing multiparameter filtrations are fundamentally different than ours.In [12], point clouds were extracted from the grayscale images, and multiparameter filtration was constructed by distances and densities among those points.Our proposed method uses grayscale images and their smoothed images directly.There also have been studies that use one-parameter persistence to study other types of biological tissues, such as self-organization of biological tissues [73], fluorescence microscopy images [74], skin [75,76], tumors [77], and other types of data [78].
There have been other mathematical models to analyze mitochondrial networks.For example, MiNA [79] is a widely used image macro tool to analyze mitochondrial networks and offers nine features to quantify mitochondrial networks.However, many of these approaches required manual thresholding.Different thresholds may result in different outcomes, and the choice of threshold is often subjective.We have addressed this limitation with our MF-PH analysis technique, which takes all threshold levels into consideration and is robust to noise.In particular, the proposed metrics-normalized Betti curves, the size distribution of Betti numbers, and the connectivity index-are scale-independent and dimensionless.Our approach is different than other groups who have developed automatic methods to segment mitochondria by a convolutional neural network [80][81][82][83].
From a clinical and scientific standpoint in bioimaging, obtaining microscopy images poses a significant challenge due to their scarcity and the complexities associated with extended experimental durations, labor-intensive processes, and resource-intensive requirements.Furthermore, images from distinct experiments often hinge on varied experimental assumptions and measurement outcomes, constraining the availability of consistent training samples and impeding the efficacy of machine learning mechanisms.
In light of these specific challenges, the employment of MF-PH methods emerges as a valuable approach.These methods offer a geometric and topological elucidation of mitochondrial networks, providing a more intuitive explanation and comprehensible understanding within the realm of medical theory.One limitation of the present study was that the proposed MF-PH framework was used on a limited sample size for each condition.Instead of cell lines, we utilized primary cells directly purified from our optineurin transgenic mice.Cells were purified from the embryos of transgenic mice.Due to the limited yield of transgenic knockout mice and the high cost of breeding these mice, our sample size was limited but sufficient to demonstrate the utility of our approach in this preliminary work.An anticipated future direction involves the utilization of our proposed MF-PH framework on an expanded set of microscopy images derived from diverse experimental conditions.
Our mathematical model and approach toward analyzing a complex biological structure using persistent homology will be useful in biomedical research and therapeutic development.In our paper, we have provided an example by applying this approach toward studying mitochondrial networks.Furthermore, we have created a "connectivity index" that summarizes all of these changes to allow quantitative comparison.Mitochondria structures take on complex shapes ranging from round or ovoid structures to long tubular networks.They can be individually isolated or connected through tubules in a branching pattern.This complex morphology presents a challenge for conventional image analysis methodology and software.However, using our approach, changes in the mitochondrial network in response to genetic mutations, aging, environmental insults, disease pathology, or drug treatment can be quantitated and compared.Beyond the mitochondrial network, further work will be required to apply our approach toward studying other complex biological structures such as endoplasmic reticulum, nuclear and cellular membranes, neuronal connectivity, and complex molecular structures.

Fig 1 .
Fig 1. Microscopy images of mitochondrial cells.Panel A: Microscopy images of mitochondrial cells in classes WT and KO, along with images of cells treated with oligomycin and antimycin A (OA).Panel B: confocal microscopy images that will be studied in this work, where the first and second rows illustrate images of WT and KO types, respectively.For the presentation, images are padded with a black background to ensure proper alignment.The original images are included in the released materials for this work.For more details, refer to the data and code availability sections.https://doi.org/10.1371/journal.pone.0310157.g001 Fig 2 once more, Figs 1 and 3(a) show examples of grayscale images.Pixel value 0 represents black while 255 represents white; values in between 1 and 255 represent the gray color gradient from black to white.

Fig 2 .Fig 3 .
Fig 2. Illustration of a binary image and its function representation.The collection of black and white pixels for this binary image f is the 12-by-12 grid denoted by P = {(x, y) | x, y = 1, 2, 3, � � �, 12}.Values 0 and 1 represent the black and white color, respectively.The set of black pixels is denoted by K(f).https://doi.org/10.1371/journal.pone.0310157.g002 Fig 4(b), those white regions whose size are less than (4 + 1) × (4 + 1) are removed by the opening operation O S 4 ðf Þ.Similarly, in Fig 4(c), those white regions whose size are less than (10 + 1) × (10 + 1) are removed by O S 10 ðf Þ.Selecting the "correct" structural element in practice is often a challenging task performed manually by the researcher and is subjective.Instead, our proposed method utilizes information from all parameters.
Fig 5(a) as an example.Since there are three disjoint clusters of black pixels (background, two small black circles inside the white ones), β 0 = 3; since there are four clusters of white pixels that are completely enclosed by the black ones, β 1 = 4. Similarly, in Fig 5(b), we observe there are ten white regions that are completely enclosed by the black pixels, and thus, β = (1, 10).Betti numbers are known to be invariant under translation, rotation, shrinking, or enlarging (cf.[39], Ch. 20, p. 128).For instance, as shown in Fig 5(a), a white vertical bar and a white circle are different shapes and sizes, but they all count equally towards β 1 .To echo its property, consider the binary image in Fig 5(b).There are various sizes of the white regions from a size of 3 × 3 to 20 × 20, but they all count equally towards the β 1 .Lastly, we demonstrate the interplay among global thresholding, opening operations, and Betti numbers.Consider Fig 6, where g is a grayscale image of a WT cell.First, let us consider what happens to Betti numbers when we use current analysis techniques employing manual thresholding.Fig 6(b) shows the thresholded binary image at 20 and its Betti number is (334, 259).We apply the opening operation with respect to S 2 on the binary image g 20 to obtain Fig 6(c), and its Betti number is (96, 42).We observe that the binary image O S 2 ðg 20 Þ visually is cleaner than g 20 and their Betti numbers change drastically.As a comparison, Fig 6(d) shows the thresholded binary image at 80.It seems that g 80 reveals a finer structure of the cells than g 20 does.We apply the opening operation with respect to S 1 on the binary image g 80 to obtain Fig 6(e).At this point, we observe that there are two thresholding parameters in this process, and choosing different pairs of parameters results in different images and Betti numbers.

Fig 5 .
Fig 5. Sample binary images and Betti numbers of black regions.β 0 counts the number of isolated black regions (counted by yellow numbers in (a) and (b), respectively), and β 1 counts the number of white regions that are completely enclosed by the black ones (counted by red numbers in (a) and (b), respectively).https://doi.org/10.1371/journal.pone.0310157.g005

Fig 8 .Fig 9 .
Fig 8. Illustration of a bifiltration of a WT image.The bifiltration is defined as in (10), where T = 20, 40, 60, 80, 100 means the threshold of image intensity, and O means the sizes of the opening operation corresponding structuring elements S 0 , S 2 , S 4 , S 6 , S 8 .The WT image is the first in the first row of Panel B in Fig 1.https://doi.org/10.1371/journal.pone.0310157.g008

Fig 11
Fig 11 shows the normalized Betti curves in (17) for those images of WT and KO type cells as in Fig 1.Since there are 5 images (for 5 different cell samples) for each genotype, there are 5 curves in each Fig 11(a), 11(b), 11(d), 11(e), 11(g), and 11(h).We first consider the results of WT images.Fig 11(a), 11(d) and 11(g) show the normalized Betti curve at the 1-st, 4-th, and 6-th row of diagram (10): r H;0 1 , r H;3 1 , and r H;5 1 , respectively.Recall that each row of the diagram (10) indicates a different level of opening/smoothing where the 1-st row indicates no smoothing, the 4-th row indicates smoothing with O S 3 , and the 6-th row indicates smoothing with O S 5 .We observe the shapes of these curves are different among different WT cell samples.This suggests that mitochondrial morphologies in WT cells exhibit large variability.Second, we consider the results of KO images.Fig 11(b), 11(e), and 11(h) show the normalized Betti curve at the 0-th, 3-rd, and 5-th row of diagram (10): r H;01 , r H;3 1 , and r H;5 1 , respectively.We observe that in (b), the shapes of these curves are similar among different sample images.

Fig 10 .
Fig 10.Illustration of opening filtrations and the corresponding connectivity indices.To illustrate different images' fragmentation and connection structures, the binary images are derived from a mitochondrial image in class WT using morphological dilation and thresholding.CI refers to the connectivity index.https://doi.org/10.1371/journal.pone.0310157.g010

Fig 11 .
Fig 11.Comparison of normalized Betti curves for WT and KO images.The normalized Betti curves are defined as in (17), and the WT and KO images are shown in Fig 1. (a)-(b) normalized Betti curves for WT and KO respectively, r H;0 1 , along the horizontal (0-th row of diagram (10)) filtration fX t;0 g 255 t¼0 where 0 indicates no opening operation is applied; (c) Average curves over those 5 images in (a) and (b); (d)-(e) normalized Betti curves for WT and KO respectively, r H;3 1 , along the horizontal (3-rd row of diagram (10)) filtration fX t;3 g 255 t¼0 where 3 indicates opening operation with structural element S 3 is applied; (f) Average curves over those 5 images in (d) and (e); (g)-(h) normalized Betti curves for WT and KO respectively, r H;5 1 , along the horizontal (5-th row of diagram (10)) filtration fX t;5 g 255 t¼0 where 5 indicates opening operation with structural element S 5 is applied; (i) Average curves over those 5 images in (g) and (h).https://doi.org/10.1371/journal.pone.0310157.g011

Fig 12 .
Fig 12.Comparison of size distributions of Betti numbers for WT and KO images.The size distributions of Betti numbers, ϕ(D t ), are defined in (19), and the WT and KO images are shown in Fig 1. (a)-(b) ϕ(D 128 ) of WT and KO, respectively, type.(c) Average over curves in (a).−0.0980 is the rate of change from size 1 to size 4. (d) Average over curves in (d).−0.2469 is the rate of change from size 1 to size 4. https://doi.org/10.1371/journal.pone.0310157.g012 curves � � : f1; :::; 20g !½0; 1� from size 1 to size 4 which are −0.0980 and −0.2469 as displayed in Fig 12(c) and 12(d).The rate of change in mitochondria size in these curves is another way to measure differences in mitochondrial size distribution between the two cell types.

Fig 13 (
Fig 13(a) shows the boxplot of the connectivity index C t for t = 100, 101, . .., 150 for images of WT (blue) and KO (red) type.There are 5 different cell samples/images per genotype shown.We observe that the connectivity index of mitochondria WT cells is greater than the KO cells.This suggests that the WT mitochondria are more connected than the KO mitochondria.To further quantify their differences, we consider the connectivity index over all possible fluorescence signal threshold values as shown in Fig13(b) and 13(c).We observe that the connectivity

Table 1 . Notations used in this article and their descriptions.
persistence diagram of a given filtration fX i g X ðt 1 ;i 1 Þ � X ðt 2 ;i 2 Þ for any t 1 � t 2 and i 1 � i 2 .

Table 3 . Hypothesis tests of the proposed features for the sampled WT and KO mitochondria images.
[72]proposed features are demonstrated inFigs 11-13.The statistics of the hypothesis tests are the average values of the WT and KO curves over the range [a, b], where the ranges are selected with reference to Figs 11-13.The test samples consist of 5 curves each from WT and KO mitochondria images.The null hypothesis, H 0 , states that the distribution of mean values of the WT and KO curves over the range [a, b] are the same.The permutation hypothesistest is performed by the python package mlxtend[72].The p-value serves as the target value for determining the based on the sub-level set filtration on a grayscale image.This statistical test demonstrates that, without morphological operations, the persistence information of images in the WT and KO classes is statistically indistinguishable.This finding indicates that the proposed morphologybased MF can provide more detailed geometric information for analyzing mitochondrial structures.
test results, and the decision value α is set to be 0.05.The hypotheses are repeated 100 times with different random seeds, and average values and standard deviations present the results.More detailed code for the statistical experiments can be found in our GitHub repository (see the "Code availability" section).https://doi.org/10.1371/journal.pone.0310157.t003